Diagrammatic quantum field formalism for localized electrons 
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We introduce a diagrammatic quantum field formalism for the evaluation of normalized expecta- 
tion values of operators, and suitable for systems with localized electrons. It is used to develop a 
convergent series expansion for the energy in powers of overlap integrals of single-particle orbitals. 
This method gives intuitive and practical rules for writing down the expansion to arbitrary order of 
overlap, and can be applied to any spin configuration and to any dimension. Its applicability for sys- 
tems with well localized electrons has been illustrated with examples, including the two-dimensional 
Wigner crystal and spin-singlets in the low-density electron gas. 
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I. INTRODUCTION 

In the last two decades considerable effort in the theory 
of electronic structure has been focused on the develop- 
ment of methods where the time for computing ground 
state properties scales linearly with the size of the sys- 
tem, referred to as 0{N) methods, N being the num- 
ber of electrons in the system^i A standard approach 
there is to make use of localized one-particle electron or- 
bitals, and to circumvent their orthogonalization though 
various strategies and approximations in the subsequent 
energy minimization. Indeed, orthogonalization involves 
computationally intense algorithms; it is particularly im- 
practical for geometrical optimizations, and it is actually 
intractable when the iV-body electron wave function is 
to be written as a linear combination of Slater determi- 
nants made from different single-particle orbitals (i.e. for 
a general spin state). 

The use of non-orthogonal orbitals, on the other hand, 
poses its own difficulties, because the antisymmctrization 
of the many-body wave function in this case introduces 
terms the magnitude of which increases as iV, N"^, N^, etc 
(leading to the well known orthogonality catastrophe). 
Thus, expectation values of operators can diverge in the 
thermodynamic limit, N — > oo. It is often the case, in 
particular for an arbitrary spin state, that there is no 
transparent or/and systematic way of dealing with such 
problems. As a result, the approximation eventually used 
may violate, for example, even the charge neutrality of 
the system and thus lead to errors that also increase with 
its size. 

In this paper, we develop a diagrammatic formalism 
to deal with such problems which can be applied for any 
spin configuration and in any dimension. We use it to 
derive a linked cluster theorem for the evaluation of ex- 
pectation values (the energy is discussed in particular) in 
terms of a convergent series expansion of overlap integrals 
of single-particle orbitals. The diagrammatic language 
is introduced by direct analogy with that of standard 
field theory. The parallel is indeed interesting, bearing 
in mind that the case of strongly localized electrons con- 
sidered here is the opposite limit of spatially uniform sys- 



tems, the traditional domain of many-body perturbation 
theory. The equivalent of the Feynman propagator will 
be seen to be the overlap integral, S, the single particle 
orbitals correspond to vertices in the diagrams, and an 
n-body operator introduces n external points. All dia- 
grams are then calculated in terms of closed loops con- 
necting the external points. Despite these similarities in 
language, however, the linked cluster expansion and the 
resulting diagrammatic rules here are quite different from 
those in standard field theory. 

Consider now a neutral system consisting of Ng elec- 
trons and Ni ions (or a uniform, positive adn rigid back- 
ground, in the case of a jellium model) in a volume V. 
The Hamiltonian of this system is given in atomic units 
by 



H = 



N„ -a m p2 

Ell + J 
^ 2M, 



1 fdr /dr'^- 

2 ./ ./ r - r' 



p(2)(r,r') 



-l-pf (r,r')-2pW(r)pf)(r') 



(1) 



where the indices e and i refer to electrons and ions, 
respectively, pj is the momentum operator of electron j, 
Pj and Mj are the momentum operator and the mass 
of ion j, and p^^-'and p^^^are the one- and two-particle 
density operators defined respectively by 



pW(r) = 5]5(r-r,), 



(2) 



and 



=i(2)(r,r') = p^^^r)p^^\r') - b{v - v')p^'-\x'). (3) 



A standard approach to solving the eigenvalue problem 
for this system is, as a first step, to find the solutions of 
the electronic problem in the clamped nuclei approxima- 
tion, where the ionic momenta are set to zero and their 
coordinates frozen. The Hamiltonian for this problem is, 
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from (P), 



^2 2j J |r-r'| 



[/5(2)(r,r')-2p«(r)pfc(r')] +C/6, (4) 



where pb (r) is the classical density of the positive charge 
(ionic or that of a uniform rigid background), Ub is its 
self-energy, and we have simplified the notation by drop- 
ping the subscript e from quantities referring to the elec- 
trons. 

In what follows we discuss the evaluation of the ground 
state properties of a system described by the Hamiltonian 
(j4]), and more specifically, the quantity 



E 



(*l*> 



(5) 



with l^f) an TV-electron trial state constructed from local- 
ized single-particle spatial orbitals centered at positions 
{Ri}. Though U) and (|4]) are formally independent of 
spin, we will nevertheless also allow for an arbitrary spin 
configuration (i.e. correlation and even order) which can 
be specified by an appropriate set of spin orbitals. 

Localized orbitals here mean that they diminish rapidly 
away from the localization centers {Ri}. The limit where 
space can be divided into regions each occupied only by 
a single one-particle function corresponds to the semi- 
classical limit where the spin configuration becomes irrel- 
evant and l^*) is a product of single-particle states. When 
this is not the case, antisymmetrization of the many-body 
wave function and the resulting exchange effects become 
an important issue in determining the structural phase 
of the ground state. 

Moving away from the semi-classical limit, and when 
the space orbitals are not orthogonal, requires a neces- 
sity to introduce terms in both the numerator and de- 
nominator of ([5]) that go as ^ 0(S'^"A^"), where S is 
an overlap integral between one-electron wave functions, 
and n = 0, 2, 3, . . .. The resulting series are obviously 
divergent as iV — > oo irrespective of how small but finite 
S is. Here, we will deal with these problems by viewing 
the overlap effects as a formal "quantum perturbation," 
which introduces some scattering of the single-particle 
amplitudes. The normalization of ([5]) is then achieved in 
a diagrammatic approach without an explicit inversion of 
an overlap matrix or a requirement to introduce a cut-off 
radius for the locahzed functions. The topology of the 
connected diagrams that give a convergent and finite ex- 
pansion for the energy (per electron) will be determined 
by the set {R^}. 

The remainder of the paper is organized as follows: 
In Section |TT] we summarize a quantum field theoreti- 
cal notation, used previously by van Dijk and Vertogen^ 
and later by Moulopoulos and Ashcroft^ for describing 
Wigner crystals. All matrix elements relevant for com- 
puting the energy are constructed from products of field 



operators. Their anticommutation relations are then 
used in Section |lIT] to develop a diagrammatic language 
for evaluating the matrix elements. In Section IIVI we 
show that the taking of a ratio of matrix elements leads 
to a linked cluster expansion. First, an algebraic expan- 
sion is obtained by generalizing a mathematical device 
used by Abarenkov^ in the context of a valence-bond 
method. Next, the new formalism is used to prove rig- 
orously that the expansion is convergent and is topo- 
logically equivalent to linked clusters of closed-loop di- 
agrams. A recipe and an example for calculating the 
energy are presented in Section [V] Further applications 
and uses of the method are discussed in Section IVll 



II. QUANTUM FIELD THEORETICAL 
NOTATION 



In the formalism of second quantization (requiring 
specification of an initiating set of single-particle states) , 
the kinetic energy and the density operators in ^ can be 
written in the forms (atomic units are used throughout): 



T : 



E_ t 

k.s 



pW(r)=^V'l(r)^.(r), 



(6) 



(7) 



and 



p(2)(r,r') = Y,4ir)4,ir')^Ar')Mr), (8) 



where cj^. ^ and Ck,s are respectively creation and anni- 
hilation operators for an electron in a plane wave state 
with a wave vector k and spin s, and and V's(r) are 

the usual field operators, i.e. 



k 



(9) 



which create and annihilate a Fermion with spin s at 
position r. A general state of the system assumes the 
form 

\^) = ^ dri---drN F{ri,si;...]rN,SN) 



xV4(ri)---C(r^)|0)- 



(10) 



Here |0) denotes the vacuum state, and the antisym- 
metrization of the wave function is implicitly built into 
(fTU)) through the anticommutation relations of the field 
operators, namely 



(11) 



and 



{V'.(r),V.'(r')} = {V'l(r),VUr')} = 0. (12) 
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In variational terms the problem is to determine the 
amplitude function F which minimizes ([S]). Because 
we want to construct the wavefunction from N single- 
particle space orbitals, the choices for F are linear com- 
binations of products of single particle functions, each 
product representing a particular fixed-spin configuration 
(the standard Hartree-Fock approximation). For exam- 
ple, the simplest ansatz, corresponding to a ferromagnetic 
(FM) state is 



N 



F(ri,si; . . . ;rAr,S7v) = H-^'^"^ 



(13) 



i=l 



where fi(ri — R^) indicates the (normalized for conve- 
nience) wavefunction of an electron localized at some po- 
sition Hi. 

Next, following van Dijk and VertogenSj we introduce 
the operators c?| and di, defined by 



7t 



dr (r)/,(r-R, 



(14) 



which create and annihilate an electron localized at po- 
sition Ri, with a one-particle function /i(r), and with 
spin Si. A state corresponding to a particular fixed-spin 
configuration, can now be written as 



1$) 



/ N 

n 

\i=l 



(15) 



and if we label all such states by, say p, a general state 
of the system can be written as a linear combination of 
terms of the form (fT5|) . i.e. 



(16) 



For example, a state corresponding to spin-singlet pairs 
of electrons will be described by;^ 



N/2 



I*) - n 



d^ 



|0). 



(17) 



Here the up and down arrows explicitly indicate the spin 
to be associated with the given operator, and it is clear 
that the electrons do not have definite spins but are nev- 
ertheless grouped in pairs where the two electrons of each 
pair always have antiparallel spins. 

From ini), (US), and it is straightforward 

to derive the following anticommutation relations for 
the newly defined creation and annihilation operators, 
namely: 



and 



{d^ 



,d,} = {dld\} 



(18) 



(19) 



where S{ij), a key quantity in what follows, is 

S{ij) = Jdrf*{r~ R,)/,(r ~ R,), (20) 

the overlap integral of two single-particle wavefunctions 
centered at Ri and R,. In addition, 



and 



{Mr),dl}^6s,sJ^{r~R^ 



{V^,(r),d.} = {Vl(r),4} = 



(21) 



(22) 



Further, if fi(k) is the Fourier transform of fi{r), then 
for a system of dimensionality D, 



(2^) 



D/2 



and 



{ck,., 4} - A^^^e-^*'-"' /,(k) (5,,,. , (23) 



(24) 



V 



{ck,.,d.} = R„4} = o. 



III. 



DIAGRAMMATIC EVALUATION OF 
MATRIX ELEMENTS 



Within the formalism of the previous section, all ma- 
trix elements of interest for the computation of the the 
energy (O assume the general form 



{0\ABC---YZ\0), 



(25) 



where A, B,C, . . . , etc. are creation and annihilation 
operators whose anticommutation relations in terms of 
localized single-particle functions have just been estab- 
lished. We now proceed to interpret these quantities as 
a sum of closed loop diagrams in a language very simi- 
lar to that of standard field-theoretical and many-body 
methods <^i5. 

We start by selecting an arbitrary labeling order of all 
distinct operators of interest; this can be done without 
loss of generality. Distinctions will be based on the label 
i for the di, r for the ^/'^(r), and k for the Ck,s operators; 
the spin label will be irrelevant. Next, we define a T- 
product of operators, T{ABC ■ ■ ■), which is a product of 
the operators A, B,C, . . ., but written in such an order 
that all annihilation operators are on the left and in de- 
scending order of their labels, all creation operators are 
on the right of the annihilation operators and in ascend- 
ing order of their labels, and the product is multiplied by 
(—1)^, where P is the number of permutations needed 
to obtain the T product from ABC ■ ■ ■ . For example. 



T{d2d3d\dldi) = {-lfd3d2did\d 



(26) 



Next, we define an A^-product (normal product) of op- 
erators, N{ABC ■ ■ ■), which is a product of the operators 
A, B,C, . . ., where all creation operators are on the left 
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of all annihilation operators, and the product is multi- 
plied by (—1)^, with P being the number of permuta- 
tions needed to obtain the A^-ordering from ABC ■ ■ ■ . 
For example, 

N{did2dl) = {-lfdldid2 = {''lfdld2di. (27) 

We can now define a pairing, or a contraction of two 
operators as 

A'^B" = T{AB)~ N{AB) 

( {A,B}, if AS is T-ordered , 
\-{A,B}, if AS is not T-ordered, ^ ^ 

and then we have the equivalent of Wick's theorem for 
our problem. This states that a T-product can be ex- 
pressed as a sum of all possible A^-products with all pos- 
sible contractions, i.e., 

T{ABC---YZ) = N{ABC ■ ■ - YZ) + NiA^B^C ■ ■ - YZ) 
+ NiA^BC ■ ■ ■ YZ) + ■■■ 
+ N{A''B''C''---Y''Z'=). (29) 

The validity of the above relation can be verified by in- 
spection, but it is also not difficult to prove by induction. 

Next, taking the vacuum expectation values of ((28)) 
and ([29]) . and using the fact that by the definition of an 
A'^-product its vacuum average is zero when the product 
contains any uncontracted operators, we have 

A'^B" = {0\T{AB)\0), (30) 

and 

{0\T{ABCD---YZ)\0) = 

{o\T{AB)\o){o\T{CD)\o) ■ ■ ■ (o|r(rz)|o) ± 
{o\T{AC)\o) {o\T{BD)\o) ■ ■ ■ (o|r(rz)|o) ± • • -(31) 

where the ± signs correspond to the parity of the permu- 
tation of the operators ABC ■ ■ ■ XYZ. As a consequence, 
any matrix element of the form (I^S)) is evaluated in com- 
plete analogy with correlation functions in field theory. 

Accordingly, we now develop a diagrammatic descrip- 
tion for such matrix elements. The operators we are deal- 
ing with have three attributes: a label associated with 
the localization center of a one-particle function (for the 
di operators these functions are the /i(r)'s, for the the 
-08 (r) operators the S{r)'s, and for the Cs,k operators, the 
/i(k)'s.); a spin orientation; and every operator is either 
of a creation or annihilation character. So, we will draw 
points to represent the set of labels of the operators (these 
points can obviously be be arranged to reflect the actual 
topology of the set {Ri}), and arrows pointing away from 
or towards them for creation or annihilation operators, 
respectively. In addition, we will indicate the spin with 
a bar across the arrows for spin up operators, resulting 
in what we will refer to as plus and minus arrows. For 



example. 



4t = '^i ^'.T = (32) 

4i = I, Al - 1, (33) 

and similarly for the V's(r) and Cs^k operators. To ex- 
tend the analogy within the language of field theory even 
further, we will call the points associated with the di op- 
erators vertices, and those associated with the ips (r) and 

k operators external points; the reason for this choice 
will become clear later. 

In this construction, a pairing of two operators is rep- 
resented by a line connecting the points associated with 
them, and having a direction determined by their order- 
ing. When the operators are T-ordered, the lines will 
follow the arrows of the points they connect. Also, be- 
cause the commutation relations of opposite spin opera- 
tors are zero, only lines connecting either plus or minus 
arrows need be considered. If it is not possible to connect 
all points in this fashion, the corresponding matrix ele- 
ment is zero. This means that if all operators are present 
in the product as creation-annihilation pairs (each point 
has two arrows, one pointing at it and one away from it), 
the resulting non-zero diagrams are a collection of closed 
loops only. 

It is now easy to see that after a couple of permutations 
the expectation values of the kinetic energy and density 
operators ([6])- ([8]) can be brought to the form (|3ip . where 
all operators are present in creation-annihilation pairs 
but all terms involving pairings between the tpsi^) and 
the Cg^k operators have canceled out. Therefore, all rel- 
evant matrix elements can indeed be evaluated as the 
sum of all possible closed-loop diagrams that can be con- 
structed by connecting all vertices and external points 
according to the rules described above. The value of each 
diagram is then a product of the values of each line con- 
necting two points, and the value of each such line is the 
anticommutation relation of the operators represented by 
the points. In addition, a sign must be associated with 
each diagram, which is given by 

I 

where is the number of external points, the product 
is over all (closed) loops in the diagram and n/ is the 
number of lines (or points) in each loop; the one-point 
loop diagrams obviously have no influence on the sign. 
The formal proof of is straightforward (e.g. by in- 
duction). 

By way of example, and to illustrate the rules de- 
rived so far, we show the diagrammatic expansion of the 
one-particle density, (\I'|/5^^^(r)|\I'), where l^") is an A^- 
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electron ferromagnetic state; thus, 



* - V - '0 

o ^ o • 

r r r m 



+ 



B 

r 

i I " 

000 

r mp 



\ I 




m 



(35) 



Here the filled points indicate summation over all ver- 
tices, and we have omitted one-point loops, which are 
equal to unity. In the semi-classical limit, or if the single- 
particle functions are orthogonal, only the first diagram 
above remains. The overlap-order of a diagram increases, 
and their values diminish exponentially with increasing 
size of the loops. The presence of disconnected loops is 
generally what causes such quantities to diverge in the 
thermodynamic limit. As with other many-body meth- 
ods, this problem is removed by the normalization of the 
expectation values, which leads to the equivalent of a 
linked cluster expansion. 



IV. CONSTRUCTION OF A LINKED CLUSTER 
EXPANSION 

Let I A'') be a product (or a linear combination of prod- 
ucts) of iV creation field operators acting on the vacuum 
state. As we showed above, can be thought of 

as a sum of all possible closed-loop diagrams connecting 
some representative N points. Then, we can write: 



(36) 



etc., or more generally, for the m"^ derivative of Qjssity. 
g(^)(0) = m! C{N-m,n2,...,nN), (42) 



subject to the constraint ^ knk = m. 

k=2 

Furthermore, we can also define a function, associated 
with Qwit), by 



Rn{t) 



QN+n{t) 
QN{t) 



(43) 



and think of the original N points as vertices represent- 
ing one-particle functions, and the additional n points 
as external and representing an n-body operator O^"^. 
Then, within this construction, 



(*|6(")|*) 



= Rn{l) 



diagrams with lines 
connecting external points 



(44) 



and a diagrammatic expansion of the above can be 
obtained by considering the Taylor expansion of R{t) 
around t = 0, namely, 



oo rn -. 



m=0 i=0 



1 



(m-i)! \QN+n{0) 



(m—i) 



(45) 

For our further discussion it will be convenient to denote 
by Vi the value of all diagrams that can be constructed 
out of ony i vertices, and that do not contain 1-point 
loops. Similarly, Xi will indicate all such diagrams, but 
where the i points may include external points. Clearly, 



y, = ig(;)(0), 



(46) 



and 



where C(ni, . . . , njv) is the class of all diagrams contain- 
ing exactly n\ 1-point loops, n2 2-point loops, and so on. 
Next, let us define a generating function, 

QN{t)= Y C(ni,...,n^)t^-"\ (37) 



of a standard continuous variable, t, and because 

N 

J2 knk = N, it is clear that 

fe=i 

Qn{0) = C(Af,0,...,0), (38) 

Q'n{0) = 0, (39) 

g^(0) = 2!C(A^-2,1,0,...,0), (40) 

Q%{0) = 3!C(Ar-3,0,l,...,0), (41) 



A,: = - 



.(0). 



(47) 



A. Fixed-spin configuration 

First, we consider the case when |A^) is a single prod- 
uct of dj operators. In such a fixed-spin configuration, 

(\ (jn—i) 
Q^+„(0) J can be decom- 
posed simply as 



1 



1 



(m—i) 



(m-i)! \QN+n{0), 



(48) 



31- --Jn 
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where ji H hjn 



m — i. Then, because i?o(l) = 1; 



we have ^ XiWm~i — 0, and it foUows that: 

1=0 



^"(1) = ^ ^ - 

m—0 i—0 
= 1 + {V2 - X2) + {V3 - X3) 

+ [V4 -Xi- 6X2(^2 - ^2)] + • • • (49) 

This expression now represents a convergent hnked clus- 
ter expansion. It is easy to see that the second and third 
terms above are simply two-point and three-point loops 
involving the external points. The higher order terms 
are more complicated, but they are equivalent to chained 
loops connected to the external points. For example, the 
fourth term is a sum of all four-point loops involving ex- 
ternal points and a product of two two-point loops, either 
chained or not, but connected to two different external 
points. To prove that (|49)) is indeed convergent and to 
illustrate the diagrammatic rules, we present a second 
construction for the linked cluster expansion. 

Thus, we will now consider only the case when n = 1 
and introduce the following notation: Lri-^...i,r, will denote 
a single loop connecting all points labeled r, ii, . . . 
\N — {ii . . .im}) will be a state obtained by removing 
operators d]^ . . . d -^^ from |iV), Ri-^^...i^ will then be de- 
fined by the ratio 



Rj 



and Dj 



(jV^{n...z„}|jV-{ii...i„J) 

{N\N) 



(50) 



will denote the subset of all diagrams from 
{N\N), which have at least one non-unity loop connected 
to any of the points {ii . . .im}, divided by {N\N). Then, 
i?i(l) can be expanded in the following way: 

1 ^ 

^l(l) = Lr -\ T / , Lri^,,,i Ri i (51) 

to! ^ — ' 
m—\ 

where summation over repeated indexes is implied, and 
l/m! is to take account of repetitions. 

To every Ri^,,,i^ term, we now add and subtract 
Di-i^...i^ leading to 



i?l(l) — Lr H : Lri^,,,i 

to! ^ — ' 

7n—l 

1 ^ 
m ' ^ — ^ 



(52) 



Subsequently, the terms Di^,,,i^^ can be decomposed into 
products of loops connected to the points {ii . . .im} and 
ratios Ri^...i^.... By repeating this procedure, we are 
building a chained structure of loops connected to the 
external point, r. Every repetition contributes a minus 
sign and exactly one more (surviving) loop to the chain. 



To see the latter, consider a particular element, say 

N-l 

N-2 



The first two terms above add one, while the third adds 
two loops to the cluster at ii and «2- However, if the 
recursion procedure is applied to the second term once 
again, it will lead to two sums equal to the third term in 
([55)) , but with opposite sign. Thus, each step of the ex- 
pansion contributes exactly one loop to the linked cluster 
that survives subsequent iterations, a minus sign to the 
diagram, and importantly, increases its order by S^. 

The expansion operation, (j52p . can be applied to all 
members of (jSip any number of times, M, for any given 
N, until we obtain a sum of all possible loops involving 
the external point, r, chained to them 0, 1,2,3,... con- 
nected loops (with repetitions) involving the N vertices, 
and a remaining leading term of the order 



0(5 



2{N+M) 



1 



{N\N) 



(54) 



The construction can be generalized to the case with rig 
external points by noting that there will be then sim- 
ply He such linked clusters connected to the external 
points, or alternatively it can be seen by considering that 
i?„(l)=i??(l), e.g. 



{N + 2\N + 2) _ {N + 2\N + 2) {N+ 1\N + 1) 



{N\N) 



(A^ + llA^+l) {N\N) 



(55) 



Summarizing to this point, the normalized expecta- 
tion value of an n-body operator equals the sum of all 
diagrams where n = Ue external points are connected 
by a single loop to linked clusters of loops connecting 
vertices. Loops with lines connecting external points di- 
rectly are not permitted, while any powers of vertex-only 
loops are allowed. We already saw that an nj-point loop 
picks up a sign (—1)"'"^. In addition, the construction 
of the linked-cluster expansion shows that the addition 
of every new loop alternates the sign, so a diagram with I 
loops, of which le connect external points, and Ue exter- 
nal points has to be multiplied also by (— 1)"'=(— 1)'^'=. 
Altogether the result is that a diagram with a total of Ni 
lines connecting distinct points has a sign given by 



-1 



(56) 



Continuing with the example of the one-particle den- 
sity in a single-determinant many-body state, the nor- 
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malization of (|35[) now gives: 



r 


r 


j 

A,- 


j 

-A 



A 



/(J) 



(57) 



Here, as before, filled dots indicate summation over all 
vertices, different labels mark distinct points, and the 
label in the last two diagrams implies that the cor- 
responding point may coincide with the point j. This 
series is an expansion of the density 



where 



/' 



and for i > 1 



drpo(r) =iV, 



Jdrp,{r) = 0. 



(58) 



(59) 



(60) 



Further, /9o(r) — J2i ^ ^^OP is the density in the 
semi-classical limit, or if the one-particle functions were 
orthogonal; pi{r) = 0; and p2(r), P3(r), P4(r), . . . are the 
terms in curly brackets in (j57p , every one of which rep- 
resents a different order of overlap. It is easy to see that 
they indeed satisfy (|60|) . because 



' o o J 



dr 



= I , 



J ■ 



(61) 



loops. With the diagrammatic formalism, on the other 
hand, it is intuitively straightforward to maintain charge 
neutrality by grouping all diagrams involving a given set 
of vertices. 

As a second example, the diagrammatic expansion of 
the two-particle density, p(^^(r,r'), consists of the prod- 
uct (r)p(^) (r'), which can be obtained from ([57)) . and 
supplemented by another part with diagrams where a sin- 
gle loop is associated with both external points, r and r'. 
Both parts contain overlap-dependent diagrams giving 
rise to exchange-correlation effects. Those coming from 
(i-)p(i) (i-') arise solely from the non-orthogonality of 
the one-particle functions; they are sometimes called in- 
direct exchange terms and are usually responsible for the 
molecular bonding (not in a ferromagnetic, but spin- 
paired state, of course). The diagrams where both ex- 
ternal points are linked with a single loop give rise to the 
so called direct exchange, and some of them, for example 



(-) /'i. 

.'a' 
J 



(62) 



survive even if the one-particle functions are orthogonal. 
Notice that direct exchange comes from parallel spin cor- 
relations, and indeed, we cannot form a loop such as (j62p 
(even with more vertices) so long as any two electrons in 
it are in an antiferromagnetic arrangement. This is not 
the case for the indirect exchange, where r and r' are in 
separate loops. The qualitative differences between the 
direct and indirect exchange can also be seen from the 
fact that same-order overlap diagrams representing the 
two terms have opposite sign (see ([55|). For instance, 
compare with 



(63) 



Thus, the diagrammatic language accurately captures the 
well known fact that the ground state electronic structure 
is often determined by the competition of the two types 
of exchange. 




Each of the diagrams in the curly brackets in (j57p repre- 
sents a localized effective exchange charge and they can 
be grouped in pairs forming electric dipoles. The terms 
Pi(j), for i > 1, actually include summation over all ver- 
tices, and therefore represent higher order multipoles, e.g. 
P2(r) is a quadrupole. 

If we now return to the expansion in (|45p and com- 
pare it with ([FT)) , we see that po(r) is given by X2 — ¥2- 
However, the term — V3 gives only the single 3-point 
loop (with a minus sign) in P2(r), so if (|45p is truncated 
at this point, charge neutrality in the system will be vio- 
lated. The required neutralizing part in ^2(1") comes from 
the next term in ([45]), which contains products of 2-point 



B. General-spin configuration 

Dealing with a general-spin configuration means con- 
fronting the fact that |A^) must be a linear combination 
of state vectors, each one written as a product of N op- 
erators. 



M 



\N)=J2b,\N) 



(64) 



here the &p's being arbitrary constants. Therefore, we 
have to consider different configurations resulting 
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from {N\N), each of them with N points but with dif- 
ferent sets of arrows, representing creation and annihi- 
lation operators and their spins. The difficulties that 
now arise are related first, the fact that a given dia- 
gram can be present in more than one configuration, 
and second, to the consideration that not all diagrams 
can be constructed in all configurations. For example, 
the C{N, 0, . . . , 0) class diagrams, which are simply (and 
only) one-point loops, exist only in the p{N\N)p config- 
urations, and there are M of them. Then, 



(65) 



rather than unity, so here we are obliged to keep a tally 
even of the one-point loops. 

As a result of all this, formulating the expansion rules 
by following the linked-cluster construction outlined be- 
tween Eqs. ([50]) and ([5i)) might appear to become quite 
cumbersome for a general state because of the required 
book-keeping, even though there are no qualitative dif- 
ferences with the single-determinant case. However, be- 
cause the ferromagnetic state leads to a complete set of 
diagrams for a given set of points {R;}, we can reason- 
ably expect that the linked-cluster expansion for a gen- 
eral state can be obtained from that of a ferromagnetic 
state (e.g. (|57p') by multiplying every term in it by a co- 
efficient related to the frequency of occurrence of its ele- 
ments over all spin configurations resulting from (A^|iV). 

This conclusion can be verified by examining Eq. (|49)) . 
ft is still valid when |A'') is in the general form how- 
ever, the Wm-i's, which were previously given by (|i5)) . 
now contain QN^nifi) to the power m — i + I in the de- 
nominator, namely: 



Wm-i 



1 



E 



X 



31 



X 



3m-i 



QjV+ri(0) 



, (66) 



where QAf+„(0) is also given by The meaning of 

the ViS and XiS also changes; while in Section IfV Al 
they were equal to the two-or-more-point loops that can 
be constructed out of i points, now we have to sum over 
all configurations coming from {N\N) where these same 
loops can be formed and where the remaining N —i points 
form one-point loops (i.e. they represent a fixed spin 
state, {N — {i}\N — {«})). In practice, the latter con- 
dition actually greatly simplifies the calculations, as will 
be demonstrated in an example below. The final result 
therefore is that the expansion (|49)) . and consequently 
the diagrammatic rules derived in Section IIV Al remain 
the same for the general case, but now every loop carries 
a coefficient, equal to 



{pp'} 

M 



(67) 



where the sum over {pp'} is over all configurations 
p{N\N)pi where (1) the given loop can be formed, and 



(2) all remaining points are of definite spins, i.e. either 

or . These coefficients can be thought of as 

weights of the various loops, and in the case when all the 
6p's are equal to unity, they are simply the fraction of all 
configurations in which the given loop diagram appears. 

Continuing with the example of the one-particle den- 
sity, the generalization of l|57p is now: 



r 




Cijl 



A, - ,A, 



o 
r 



+ CijCii < j 



J 




(68) 



Here, the coefficient associated with the two-point loop 
connecting r and i is unity for normalization reasons, and 
the terms in brackets must have the same coefficients in 
order to preserve the charge neutrality; both of these 
statements are actually easy to verify explicitly. So, to 
obtain the expansion to the given order of overlap for a 
particular spin state, it is only necessary to determine 
the coefficients for two- and three-point vertex loops, Cy 
and Cijl respectively. We will now show how this is done 
with the example of a spin-singlet paired state (fT7|) . 

The wave-function in (jl7p is a linear combination of 
M = 2^/^ products of field operators with bp = ±1, so 
Ep bp = 2^/^. If we pick a particular spin pair, it leads 
to four types of configurations in (5'|5'): 



(i) 



(ii) V 

(iii) - \/ 

(iv) - 

To determine Cy , we have to count all configurations 
where we can form a two-point loop out of i and j, and 
form one-point loops of the remaining points. Thus, if 
i and j belong to the same pair, they must be either in 
state (iii) or (iv), thus bringing a factor of —2. The re- 
maining N — 2 points must be either in configuration (i) 
or (ii), of which there are ^2^/^, and all of them with pos- 
itive sign. So, in this case, Cy = -2^2^/72^/2 = -1. If 
i and j belong to different pairs, both of these points have 
to be either as in (i) or (ii) . This is because the remaining 
points from each pair must form one-point loops. From 
the remaining four combinations only two survive, be- 
cause i and j must be associated with parallel spins. So, 
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the two pairs contribute two configurations, the remain- 
ing TV — 2 points, as before, give rise to ^2^/^ possible 
diagrams with only one-point loops, and we therefore find 

For Ciji, we have to consider 3-point loops; they can 
connect either 3 points all belonging to different pairs, or 
3 points two of which can be from the same singlet pair. 
In the former case, all points must be in configurations 
(i) or (ii), as was the case with Cij, and then they all 
have to be associated with parallel spins. There are two 
such configurations, and the remaining iV — 3 points give 
i2A'/2 j[YioTe, so the result is that c^; = j. If 2 of the 
3 points belong to the same pair, they must be either 
in (iii) or (iv). Then, in either case, the remaining point 
must be in either (i) or (ii), but not in both. So, there are 
two options each carrying a minus sign. The remaining 
A'' — 2 points give ^2^/^ eligible combinations, and in this 
case Ciji = ^. 

To summarize, we have determined that for the spin- 
singlet paired state lfT7|) . 

{—1, if i and j are in the same pair 
(69) 
i, if i and j are not in the same pair. 



and 



Ciil 



if any 2 oi i, j, I are in the same pair 
if i, J, I are from different pairs. 



n n 

Result ([69)) is in agreement with Ref'sy and|j, but here 
it is obtained in a quite different way; and with (j70p we 
are going one step further, as we already have the next 
term in (j68p without further effort. In fact, the expansion 
([68)1 has seven terms (if we open the brackets), however, 
with the diagrammatic language it is easy to see first, 
that only two of their coefficients are unique, and next 
to determine them. 



2. Determine the coefficients associated with the spin 
configuration for all diagrams up to the required 
order of overlap. The order of overlap of a diagram 
is usually equal to (but may be higher than) the 
number of interconnected vertices in it. 

3. Form all connected, topologically non-equivalent 
and non-zero diagrams with one and two external 
points up to the required order of overlap following 
the rules described in Section HVl 



4. Determine the signs and symmetry factors (multi- 
plicity) of all diagrams. 

5. Group diagrams involving the same vertices; each 
group represents either a direct or an indirect (with 
zero net charge) exchange term. 

6. With each solid line associate an overlap integral: 
i . . j -> S{ij) = JdvMr - R,)/, (r - R,), 

generally assumed to be small. 

7. With each external fc-point associate a kinetic en- 
ergy term: 



•J 



T{^J) = y'dk^/,(k)/,(k)e-*(^-^^) 



i. With each pair of external points r and r' associate 
an interaction energy term: 



■ • • • 



U{ij, kl) 



V. ENERGY CALCULATION 

In this section, we demonstrate the use of the diagram- 
matic technique for evaluating the energy of a system 
with localized electrons. First, we formulate general rules 
for such calculations, and then we apply them to a prac- 
tical example. 



A. Diagrammatic rules 

To calculate the energy with the help of the diagram- 
matic language, we adhere to the following procedure: 

1. Specify the localization points, {Ri}, for the single- 
particle functions, fi(j), and decide the required 
order of overlap. 



where here /i(r) stands for /^(r — R)), etc. 
9. Sum over all vertex points, and m. 

The advantage in following this procedure is that 
all exchange-correlation terms originating from the non- 
orthogonality of the one-particle orbitals can be easily 
pre-summed, thus reducing the complexity of the prob- 
lem to that of one with orthogonal orbitals. The com- 
putational cost is then limited by the efficiency for the 
evaluation of the Coulomb repulsion integrals U{ij,kl). 
Their computation can be carried out with existing algo- 
rithms that scale linearly beyond a given N, for example, 
the linear scaling methods developed by Schwegler and 
Challacombei^, for computation of the U {ij, kl) integrals 
based on multipole expansions. 
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B. Example: The two-dimensional Wigner crystal 

As an illustrative example of the application of the 
procedure described above, we consider the case of the 
ground state of a 2-D Wigner crystal^ (WC) where the 
electrons are localized on a hexagonal lattice in the pres- 
ence of a uniform, rigid and neutralizing background. For 
A'' electrons in an area A, the background charge density 
is pb = N/A = l/7rr^, which also defines the dinien- 
sionless density parameter r^. Quantum Monte Carlo 
calculations have predicted that the 2-D WC exists for 
rs > 37= The hexagonal lattice has primitive vectors 



ai = a(l, 0), and 



a2 = ^(l,V3), 



(71) 



where a = y "^^s is the lattice parameter, and the elec- 
trons are localized on lattice sites = n^ai -f miSL2. 
With each electron, we associate a normalized Gaussian 
(trial) wavefunction in 2-D with width a, i.e. 



/(r-R.) 



J_g-(r-R.)V2T^ 

/ 9 



(72) 



The choice of Gaussians here is justified because the po- 
tential around the equilibrium positions of the electrons 
is close to harmonic.^" Finally, we will restrict our dis- 
cussion to the antifcrromagnetic (AFM) state - a spin- 
frustrated structure with alternating lines of up and down 
spins, e.g. an electron localized at Ri — niSii + niia.2 will 
have a positive (negative) spin if mi is even (odd). 

With these preliminaries, we can now proceed to cal- 
culate the energy per electron. The overlap integral be- 
tween one-particle functions centered at R^ and Ro is 



(73) 



where Ry = Ri — Rj. Typical values for a can be 
estimatedii to be less than a/4 and, therefore, the near- 
est neighbor (n.n.) overlap is S" = 'S'(a) < e^^ w 0.37. 
Since S"^ « 0.02 and ~ 0.007, inclusion of diagrams 
up to 0{S^) will guarantee a better than 1% precision 
in the calculation of the total energy. The next n.n. 
(n.n.n.) distance in the triangular lattice is v^a, which 
means that the n.n.n. overlap integral is S{'\/ia) = . 
Therefore, for the required precision we need to consider 
only n.n. overlaps because the two-vertex diagrams are 
of order ^^(-v/Sa) — and the three- vertex are of order 
S{a)S{'\/ia)S{a) — 5 when they involve n.n.n. overlaps. 

The relevant coefficients associated with the spin con- 
figuration are 



and 



Ciil 



^Si ,Sj ^Si .Si ^Sj ,Si 1 



(74) 



(75) 



where Si indicates the spin of an electron localized at 
Ri. For diagrams involving only n.n. overlaps, i.e. when 



Rij = Rii = Rji = a, we have Ciji — 0, and we can also 
set: 



1, If Ri, = ai 



(76) 



[0, ifR,, ^ai. 

Thus, the three-vertex 0(5*^) diagrams vanish, and we 
are left with only two-vertex diagrams of order S'^ and 

The relevant diagrams with one external point, to- 
gether with their signs and multiplicity factors are: 



0(^°) : 
0(5^) : 
0{S^) : 



% 

k 



x5 
k 




-3 



(77) 
(78) 
(79) 



Each of the 0(5*^) diagrams above has in principle a sym- 
metry factor of two - in the first diagram the external 
point can be connected to either i or j, and in the second 
the three-point loop can go either clockwise or counter- 
clockwise. However, this symmetry factor is taken care 
of when performing a sum over i and j and allowing rep- 
etition, e.g. = {1,2; 2,1} (but i ^ j). The 0{S^) 
diagrams have a multiplicity of three here, because in the 
triangular lattice there are three diagrams of each type, 
namely, in addition to those shown above, when I = i, 
and with a il loop instead of jl (for a different spatial 
or spin configurations we may have to write these explic- 
itly). 

Given ([77)) - (ITS)) , we associate the following kinetic en- 
ergy terms: 

T{zt) + T{ii)S\ij)~T{ij)S{i]) + 

:^[T{zt)S\i])-T{zj)S{zj)\ S\jl), 

which with the choice of Gaussian wavefunctions have a 
simple analytical form, namely 



T{zj) = S{zj)TiO) 1 - 



4o-2 



(80) 



Here T(0) = 1/2ct-^ is just the energy of a 2-D harmonic 
oscillator. Then, after summing over i, j, and we obtain 
the kinetic energy per electron as: 



T 
TV 



r(o) 



The relevant diagrams with two external points, rep- 
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resenting the electron-electron interaction energy are: 
i 



0{S°) 



1 , . , , 

2 o' o 
r r' 

^ o ' * o 



r r 



r r' 



1 '•^ 

2 cs' b 
r r' 



(82) 



(83) 



(84) 



(5 

r' 



^ cS o 
r r ' 

The factor of 1/2 conies from the symmetry with re- 
spect to exchanging r and r', as in Eq. ([4]), and it takes 
care of overcounting. Notice also that the 0{S°) term 
is the Hartree interaction, while the last terms in the 
0(5*^) and 0{S^) expansions are the direct exchange. 
The remaining terms originate from the product of one- 
particle density expansions, p^^^ (r)p*^^) (r') and represent 
multipole interactions. With the above diagrams we now 
associate matrix elements U{ij,lm), the general form of 
which can be simplified by substituting for the /(r) func- 
tions and l/|r — r'| their Fourier transforms: 

U {ij, Im) — 

^ y^ki Jdk2jdk3jdk4 /(ki)/(k2)/(k3)/(k4) 

gi(ki-Ri-|-k2-Rj-l-k3-R„+k4-R,„) 

X r. -i—. (5(ki+k2 + k3 + k4). 

|K3 -I- K4I 

Then, changing the integration variables according to: 
ki = — k, k2 = k — q, ka = — k' and k2 = k' -|- q, and 
using 



fdk/(k)/(k±q)e-*'^-^^ 

we obtain 

U {ij, ml) = 

S{ij)S{lm) U 

where here 



(85) 



R,„ 



2 r^'^ 

U{v) = d— dipe' 



TTCT^ JO 



(r^/2(T^) cos^ ip 



(86) 



is the interaction energy between two Gaussian unit 
charges with centers separated by a distance r. 



After summing ([82l) - ((84)) over i,j,l, and m, the 
electron-electron interaction energy per electron is then 
given by: 



% = i y u{R) 

+ {23^ 



(87) 



45* 



105*) ^ [C/(R) 
Rt^o 

ai 



2U[-)-UiO) 
10i7(y) -4C/(ai)-3[/(0) 



where terms involving U{a.i) and J7(ai/2) have been 
added and subtracted in order to complete the second 
sum above, and C/(0) = \/27r/(T comes from the direct 
exchange. For the total interaction energy, the electron- 
background and background-background energies have to 
be added, which together with the first term in ([57]) can 
be evaluated by the Ewald lattice summation method. 
The second sum is equivalent to the interaction energy 
of an ionic lattice with opposite charges at R and R + ai 
and again is straightforward to obtain by the Ewald con- 
struction. Finally, the remaining terms require only the 
numerical computation of two integrals such as given by 

The solution thus obtained straightforwardly here, up 
to and including fourth order in overlap, is to be com- 
pared with Refs. 0113) where similar problems are dis- 
cussed but only up to 0{S^). The procedure can easily 
be extended, if needed, to higher orders. 



VI. DISCUSSION AND FURTHER EXAMPLES 

As with other diagrammatic techniques, the benefits 
here come from a translation of an algebraic formalism 
into a more intuitive diagrammatic language. It provides 
insight helpful for dealing with various spin correlations 
and overlap effects of arbitrary order. The diagrammatic 
rules also offer guidance for calculating normalized ma- 
trix elements in a most efficient way for a desired accu- 
racy, while also preserving charge neutrality in the sys- 
tem. Violation of charge neutrality as a result of an ap- 
proximate treatment of overlap effects may become a se- 
rious issue depending on the system size. We will illus- 
trate this problem with an example of a system of spin 
singlets, which has been studied previously in the context 
of a low-density electron g cLS ^^^^ but also has relevance 
for molecular systems. 

Consider therefore a collection of n spin-singlet pairs of 
electrons {N — 2n) . We will assume that the separations 
between pairs are sufficient so that interpair overlaps can 
be ignored. Without loss of generality we will also set 
all intrapair separations to be the same and equal to a 
(the relevant overlap integral being S = 5(a)). For sim- 
plicity, we will examine only the exchange corrections to 
the kinetic energy; these can be easily determined to all 
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orders of intrapair overlap. With this construction, there 
are three types of diagrams relevant for the kinetic en- 
ergy expansion. They are given by (j77p and (j78p . with 
the only difference being that the sign of the 0(5*^) di- 
agrams must be changed, since the corresponding coeffi- 
cients Cij , as given by (|76p. are equal to —1. To obtain 
the expansion to all orders of S, we have to multiply ((75|) 
by closed ij loops of all powers resulting in a geometric 
progression. The kinetic energy per electron is thus given 
by: 



^ = [T{Q) + ST{a)]{l-S' + S' 



r(0) + ST{a) 



1 + 5*2 



(88) 



If on the other hand we decide to first truncate the 
expansions of (^'|T|\1/) and (^'l^') to a particular order, 
say 5*2, and then compute T/N, the result is: 



T 
N 



T{0) + in-l)S^T{0) + ST{a) 



The difference between 



1 + 



and 



ST 
'N 



T{0)-Tia)/S 
l + {n + l)S^ + nS^ 



(89) 



(90) 



which shows that if {n — 1)S'^ ~ 1, the error in 
is comparable to the leading order exchange term. In 
fact, if 71 oo, jlSl) gives T = iVr(O), i.e. 100% error 
in the exchange energy. Even if the analysis is carried 
for a central pair and only its nearest neighbors, in a 
typical crystalline arrangement n ^ 10, so that there is a 
very stringent limitation on the allowed overlap, namely 
52<l/n-0.1. 

The diagrammatic formalism can also be used to ex- 
amine the efficiency of dealing with the non-orthogonality 
problem by introducing a cut-off radius, Rc, for the one- 
particle functions, so that /(r) = if r > Rc- If the de- 
sired accuracy is second order in overlap, we know that 
only energy terms corresponding to two- vertex diagrams 
such as 



need to be considered. However, with a cut-off approach, 
even if Rc is chosen smaller than the nearest neighbor 
distance, one end up calculating (explicitly or implicitly) 
terms of higher than the required order in overlap, cor- 
responding to diagrams such as: 

If the profile of the wavefunctions requires larger cut- 
off radius, the efficiency would diminish even further as 
terms corresponding to 3-point, 4-point,etc. loops, would 
now enter the calculations. 

As a final example of application of the diagrammatic 
formalism, we will use it to gain insight into the physics 



underlying the linear scaling density functional theory 
developed by Mauri, Galli and Car^i^iiiii^ and by Or- 
dejon et al.}^ In this approach, non-orthogonal one par- 
ticle functions are used and the inverse of the overlap 
matrix, S^-*-, entering the energy functional is replaced 
by a truncated series expansion: 



M 



S-i«Q = ^(I-S)' 



(91) 



n=0 



where I is the identity matrix, and S has components 
S{ij). In addition, the following term is added to the 
energy functional: 



drp{r) 



(92) 



where 77 is a parameter that can be freely chosen, N is the 
number of electrons in the system and p(r) is the charge 
density computed with the truncated series expansion, 
Eq. (PT|) . This method does not require explicit orthogo- 
nalization; a minimization of the energy functional of the 
non-orthogonal Kohn-Sham orbitals naturally leads also 
to orthogonalization. 

Previously, the minimization procedure has been 
shown to be convergent when the expansion (j9ip is trun- 
cated at M odd and r] is chosen to be positive. We 
can easily see the physical reason for this. In the di- 
agrammatic language, an expansion of (pT|) to odd M 
corresponds to considering only diagrams for the den- 
sity expansion where the maximum number of solid lines 
(representing overlap integrals) is also odd. The expan- 
sion of the density, given by ((57)) . shows that truncating 
the series (j9ip in this way introduces an error, which is 
equivalent to decreasing the electron charge density and 
the system becoming not neutral. The extra term added 
to the energy functional, (|92p . then represents the inter- 
action energy between a positive external field and net 
positive charge. Thus, reducing this interaction energy 
to zero, i.e. energy minimization, is only achieved when 
orthogonality is attained. 

With this physical picture in mind, it is easy to see that 
the method should also work when M is chosen to be odd 
and the parameter rj negative. Indeed, in this case the 
error introduced by the truncated expansion (j9ip leads 
to increasing, not decreasing, electron charge. But with 
77 < 0, this excess charge now interacts with a negative 
field and (|92p is again positive definite. Realizing this 
without the physical picture in mind is not straightfor- 
ward because the quantity (Q — S"-^), which is negative 
definite when M is odd (see Ref. 13; for details), is not 
positive definite when M is even. 

We note that for an infinite periodic system with a 
net charge, the long-range Coulombic potential would in 
principle lead to divergent energy. In practice, the di- 
vergence can be removed by setting the q = Fourier 
component of the interaction energy to zero - this is 
equivalent to adding a uniform potential and does not 
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lead to structural changes. The remaining part of the in- 
teraction energy coming from the artificial net charge will 
be small compared to if 77 is chosen sufficiently large, 
and will also vanish when orthogonalization is attained. 

The above discussion illustrates the utility of the dia- 
grammatic formalism to inspect charge neutrality; it is 
ensured with a proper grouping of diagrams, as shown in 
Eq. jSZl). 

VII. CONCLUSION 

We have introduced a diagrammatic formalism for the 
calculation of normalized expectation values in terms of 
convergent series expansions in powers of one-particle 
overlap integrals. It can be applied to any order of over- 
lap and for any spin configuration. The formalism has 
been introduced by analogy with conventional field theo- 
retical methods; however it is applicable for systems with 
well localized electrons. As a particular example, we have 
demonstrated energy calculations up to fourth order in 
overlap at the level of unrestricted Hartree-Fock and the 
valence-bond methods. The formalism presented here 
can give useful physical insight for the validity of other 
approaches, and potentially be used improve their effi- 
ciency. 

A possible extension of the formalism can include an 



analogy of skeletal diagrams and Dyson- like equations. 
This would be particularly useful in cases where there 
is a significant overlap among groups of electrons. In 
such cases, selected diagrams, accounting for the overlap 
among these electrons, could be summed to an infinite 
order. This possibility is demonstrated with the example 
from the preceding section, Eq. ([55|) . 

The formalism can also be readily applied for local- 
ized bosons. The only difference with the fermionic case 
is in the sign of the diagrams as expected. For bosons, 
all loops carry a positive sign as a result of the com- 
mutation relations, however, in the construction of the 
linked-cluster expansion each chained loop still brings a 
negative sign. Therefore, in this case the sign of a dia- 
gram is given by (— 1)'~'= instead of ([55]) . where I — If. 
is the number of closed loops not connected to external 
points. 
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